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Abstract 

There has been a long-standing controversy whether information in neuronal net- 
works is carried by the firing rate code or by the firing temporal code. The current 
status of the rivalry between the two codes is briefly reviewed with the recent stud- 
ies such as the brain-machine interface (BMI). Then we have proposed a generalized 
rate model based on the finite A^-unit Langevin model subjected to additive and/or 
multiplicative noises, in order to understand the firing property of a cluster contain- 
ing neurons. The stationary property of the rate model has been studied with the 
use of the Fokker-Planck equation (FPE) method. Our rate model is shown to yield 
various kinds of stationary distributions such as the interspike-interval distribution ex- 
pressed by non-Gaussians including gamma, inverse-Gaussian-like and log-normal-like 
distributions. 

The dynamical property of the generalized rate model has been studied with the 
use of the augmented moment method (AMM) which was developed by the author 
[H. Hasegawa, J. Phys. Soc. Jpn. 75 (2006) 033001]. From the macroscopic point of 
view in the AMM, the property of the A'^-unit neuron cluster is expressed in terms of 
three quantities; /i, the mean of spiking rates of R = {1/N) ^ , where denotes the 
firing rate of a neuron i in the cluster: 7, averaged fiuctuations in local variables (r;): p, 
fluctuations in global variable (R). We get equations of motions of the three quantities, 
which show p y/N for weak couplings. This implies that the population rate code 
is generally more reliable than the single-neuron rate code. Dynamical responses of 
the neuron cluster to pulse and sinusoidal inputs calculated by the AMM are in good 
agreement with those by direct simulations (DSs). 

Our rate model is extended and applied to an ensemble containing multiple neuron 
clusters. In particular, we have studied the property of a generalized Wilson-Cowan 
model for an ensemble consisting of two kinds of excitatory and inhibitory clusters. 
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1 Introduction 



Human brain contains more than 10^'^ neurons [T]. Neurons communicate information, 
emitting short voltage pulses called spikes, which propagate through axons and dendrites to 
neurons at the next stage. It has been a long-standing controversy how neurons communicate 
information by firings or spikes The issue on the neural code is whether information 

is encoded in the rate of firings of neurons {rate code) or in the more precise firing times 
(temporal code). The rate code was first recognized by Adrian 7J who noted that the neural 
firing was increased with increasing the stimulus intensity. The firing rate r(t) of a neuron 
is defined by 



where tw denotes the time window, n{t, t + tio) the number of firings between t and t + t^, 
and tk the kth firing time. It has been widely reported that firing activities of motor and 
sensory neurons vary in response to applied stimuli. In the temporal code, on the contrary, 
the temporal information like the inter-spike interval (ISI) defined by 



and its distribution tt(T), ISI histogram (ISIH), are considered to play the important role 
in the neural activity. More sophisticate methods such as the auto- and cross-correlograms 
and joint peri-stimulus time histogram (JPSTH) have been also employed for an analysis 
of correlated firings based on the temporal-code hypothesis. There have been accumulated 
experimental evidences, which seem to indicate a use of the temporal code in brain [2][8]-[TT]. 
A fiy can react to new stimuli and the change the direction of fiight within 30-40 ms [2J. In 
primary visual cortices of mouse and cat, repetitions of spikes with the millisecond accuracy 
have been observed [8]. Humans can recognize visual scenes within tens of milliseconds, 
even though recognition is involved several processing steps [9]-[Tl]. 

The other issue on the neural code is whether information is encoded in the activity of a 
single (or very few) neuron or that of large number of neurons (population code) [12] [13] . The 
classical population code is employed in a number of systems such as motor neurons [H] , the 
place cell in hippocampus [15 and neurons in middle temporal (MT) areas [16]. Neurons in 
motor cortex of the monkey, for example, encode the direction of reaching movement of the 
arm. Information on the direction of the arm movement is decoded from a set of neuronal 
firing rates by summing the preferred direction vector weighted by the firing rate of each 
neuron in the neural population |14| . It has been considered that the rate code needs the 
time window of a few hundred milliseconds to accumulate information on firing rate, and 
then its information capacity and speed are limited compared to the temporal code. It is, 
however, not true when the rate is averaged over the ensembles (population rate code). The 
population average promotes the response, suppressing the effects of neuronal variability 
[T7] [T5] . The population rate code is expected to be a useful coding method in many areas 
in the brain. Indeed, in recent brain-machine interface (BMI) or brain-computer interface 




(1) 



Tk 



tk+i — tk, 



(2) 
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(BCI) [IS][in]IlI]j the information on hand position and velocity etc. at time t is modeled 
as a weighted linear combination of neuronal firing rates collected by multi-electrodes as 
given by ^0} 

yW=5^a(w)r(^-^.). (3) 

u 

Here y(i) denotes a vector expressing position, velocity etc., r(t — u) expresses a vector of 
firing rates at time t and timelag u, and a(u) stands for a vector of weight at the timelag 
u. Equation (3) is transformed to a matrix form, from which a vector a(u) is obtained by 
the linear filter method [22] [23] with training data of y{t) and r{t — u). The predictor of a 
new, unobserved stimulus y{t) for a observed v{t — u) is then given by 

u 

Visual and sensory feedback signals like the pressure on animal's skin may be sent back 
to the brain. It has been reported that an artificial hand is successfully manipulated by 
such decoding method for observed firing-rate signals of f(t — u) [Tl][ini- A success in BMI 
strongly suggests that the population rate code is employed in sensory and motor neurons 
while it is still not clear which code is adopted in higher-level cortical neurons. 

The microscopic, conductance-based mechanism of firings of neurons is fairly well un- 
derstood. The dynamics of the membrane potential Vi of the neuron i in the neuron cluster 
is expressed by the Hodgkin-Huxley-type model given by [21] 

= Y.9n{V^.C^n)[Vr-V,)+h. (5) 
n 

Here C expresses the capacitance of a cell: V,. is the recovery voltage: gn{Vi,an) denotes 
the conductance for n ion channel (n =Na, K etc.) which depends on Vi and a„, the gate 
function for the channel n: is the input arising from couplings with other neurons and 
an external input. The dynamics of a„ is expressed by the first-order differential equation. 
Since it is difficult to solve nonlinear differential equations given by Eq. (5), the reduced, 
simplified neuron models such as the integrate-and-fire (IF), the FitzHugh-Nagumo (FN) 
and Hindmarsh-Rose (HR) models have been employed. In the simplest IF model with a 
constant conductance, Eq. (5) reduces to the linear differential equation, which requires an 
artificial reset of Vi at the threshold 6. 

It is not easy to analytically obtain the rate r{t) or ISI T(t) from spiking neuron model 
given by Eq. (5). There are two approaches in extracting the rate from spiking neuron 
model. In the first approach, Fokker-Planck equation (FPE) is applied to the IF model to 
calculate the probability p{V,t), from which the rate r{t) is obtainable. In order to avoid 
the difficulty of the range of variable: —oo < y < 6* in the original IF model, sophisticate 
quadratic and exponential IF models, in which the range of variable is — cxd < < cxd, have 
been proposed [25] . In the second approach, the rate model is derived from the conductance- 
based model with synapses by using the / — / relation between the applied dc current / and 
the frequency / of autonomous firings [26] [27] [28] . It has been shown that the conductance- 
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based model may lead to the rate model if the network state does not posses a high degree 
of synchrony 28J. 

In the rate-code hypothesis, a neuron is regarded as a black box receiving and emitting 
signals expressed by the firing rate. The dynamics of the rate ri(t) of the neuron i in a 
cluster is expressed by 



where r denotes the relaxation time, Wij the coupling between neurons i and j, and li 
an external input. The gain function K(x) is usually given by the sigmoid function or 
by the / — / relation mentioned above. One of disadvantages of the rate model is that the 
mechanism is not well biologically supported. Nevertheless, the rate model has been adopted 
for a study on many subjects of the brain. The typical rate model is the Wilson-Cowan 
(WC) model, with which the stability of a cluster consisting of excitatory and inhibitory 
neurons is investigated [injlSH]- The rate model given by Eq. (6) with K{x) = x is the 
Hopfield model [31] , which has been extensively employed for a study on the memory in the 
brain with incorporating the plasticity of synapses into Wij . 

It is well known that neurons in brains are subjected to various kinds of noises, though 
their precise origins are not well understood. The response of neurons to stimuli is modified 
by noises in various ways. Although firings of a single in vitro neuron are reported to be 
precise and reliable [32j . those of in vivo neurons are quite unreliable, which is expected to 
be due to noisy environment. The strong criticism against the temporal code is that it is 
not vulnerable to noises, while the rate code is robust against them. 

Noises may be, however, beneficial for the signal transmission in the brain against our 
wisdom. The most famous phenomenon is the stochastic resonance |33| , in which the signal- 
to-noise ratio of subthreshold signals is improved by noises. It has been shown that the noise 
is essential for the rate-code signal to propagate through multilayers described by the IF 
model: otherwise firings tend to synchronize, by which the rate-code signal is deteriorated 
[34] [35] . Recent study using HH model has shown that firing-rate signals propagate through 
the multiplayer with the synchrony [36] . 

It is theoretically supposed that there are two types of noises: additive and multiplica- 
tive noises. The magnitude of the former is independent of the state of variable while that 
of the latter depends on its state. Interesting phenomena caused by the two noises have 
been investigated [37| . It has been realized that the property of multiplicative noises is 
different from that of additive noises in some respects. (1) Multiplicative noises induce the 
phase transition, creating an ordered state, while additive noises are against the ordering 
[38] [39]. (2) Although the probability distribution in stochastic systems subjected to ad- 
ditive Gaussian noise follows the Gaussian, it is not the case for multiplicative Gaussian 
noises which generally yield non-Gaussian distribution [40]- [l^. (3) The scaling relation 
of the effective strength for additive noise given by P{N) — /3(1)/a//V is not applicable to 
that for multiplicative noise: a{N) ^ a{l)/^/N, where a{N) and (3{N) denote effective 
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strengths of multiplicative and additive noises, respectively, in the A^-unit system [46j. A 
naive approximation of the scaling relation for multiplicative noise: a{N) ~ a{l)/^/N as 
adopted in Ref . ^ , yields the result which does not agree with that of direct simulation 

(Ds) m- 

Formally, noise can be introduced to the spiking and rate models, by adding a fluctuating 
noise term on the right-hand side of Eqs. (5) and (6), respectively. ISI data obtained from 
neuronal systems have been analyzed by using various population rate methods (for a recent 
review, see Refs. [47l[48]). Experimental ISI data cannot always be described in terms of a 
single probability distribution. They are fitted by a superposition of some known probability 
densities such as the gamma, inverse-Gaussian and log-normal distributions. The gamma 
distribution with parameters A and ^ is given by 

which is derived from a simple stochastic integrate-and-fire (IF) model with additive noises 
for Poisson inputs [35], r(a;) being the gamma function. For A = 1 in Eq. (7), we get 
the exponential distribution relevant to a simple Poisson process. The inverse Gaussian 
distribution with parameters A and /i given by 



A 



A(a;-^)2 



(8) 



is obtained from a stochastic IF model in which the membrane potential is represented as a 
random walk with drift [50j . The log-normal distribution with parameters fi and a given by 



Pln{x) = -j=L=—exp 
V zvrcr^ X 



2cr2 



is adopted when the log of ISI is assumed to follow the Gaussian [5T] . Fittings of experimen- 
tal ISI data to a superposition of these probability densities have been extensively discussed 
in the literature [47]-[5T]. 

Much study has been made on the spiking neuron model for coupled ensembles with the 
use of two approaches: direct simulations (DSs) and analytical approaches such as FPE and 
moment method. DSs have been performed for large-scale networks mostly consisting of IF 
neurons. Since the time to simulate networks by conventional methods grows as iV^ with N, 
the size of the network, it is rather difficult to simulate realistic neuron clusters. In the FPE, 
dynamics of neuron ensembles is described by the population activity. Although the FPE is 
powerful method which is formally applicable to the arbitrary N , actual calculations have 
been made for N — 00 with the mean-field and diffusion approximations. Quite recently, 
the population density method has been developed as a tool modeling large-scale neuronal 
clusters [52], [53]. As a useful semi-analytical method for stochastic neural models, the 
moment method was proposed[54]. For example, when the moment method is applied to iV- 
unit FN model, original 2A^-dimensional stochastic equations are transformed to iV(2iV + 3)- 
dimensional deterministic equations. This figure becomes 230, 20300 and 2 003 000 for 
N = 10, 100 and 100, respectively. 
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In many areas of the brain, neurons are organized into groups of cells such as column 
in the visual cortex [55]. Experimental findings have shown that within small clusters 
consisting of finite number of neurons (~ 10— 1000), there exist many cells with very nearly 
identical responses to identical stimuli [55j . Analytical, statistical methods having been 
developed so far are mean-field-type theories which may deal with infinite-size systems, 
but not with finite-size ones. Based on a macroscopic point of view, Hasegawa 56 has 
proposed the augmented moment method (AMM), which emphasizes not the property of 
individual neurons but rather that of ensemble neurons. In the AMM, the state of finite 
iV-unit stochastic ensembles is described by a fairly small number of variables: averages and 
fluctuations of local and global variables. For A''-unit FN neuron ensembles, for example, 
the number of deterministic equation in the AMM becomes eight independent of N. This 
figure of eight is much smaller that those in the moment method mentioned above. The 
AMM has been successfully applied to a study on the dynamics of the Langevin model and 
stochastic spiking models such as FN and HH models, with global, local or small-world 
couplings (with and without transmission delays) [57j-|61|. 

The AMM in Ref. [56| was originally developed by expanding variables around their 
stable mean values in order to obtain the second-order moments both for local and global 
variables in stochastic systems. In recent papers [46] [62], we have reformulated the AMM 
with the use of FPE to discuss stochastic systems subjected to multiplicative noise: the 
FPE is adopted to avoid the difficulty due to the Ito versus Stratonovich calculus inherent 
to multiplicative noise |63j . 

An example of calculations with the use of the AMM for a FN neuron cluster subjected 
to additive and multiplicative noises is presented in Figs. l(a)-(e) [64]. When input pulses 
shown in Fig. 1(a) are applied to the 10- unit FN neuron cluster, the membrane potential 
Vi{t) of a given neuron depicted in Fig. 1(b) is obtained by DS with a single trial. It has the 
much irregularity because of added noises. When we get the ensemble-averaged potential 
given by V{t) ~ (1/A^) X^i ^i(^)> the irregularity is reduced as shown in Fig. 1(c). This is 
one of the advantages of the population code [17]. The results shown in Figs. 1(b) and 1(c) 
are obtained by a single trial. When we have repeated DS and taken the average over 100 
trials, the irregularity in V(t) is furthermore reduced as shown in Fig. 1(d). The result of 
the AMM, /i(t), plotted in Fig. 1(c) is in good agreement with that of DS in Fig. 1(d). 

The purpose of the present paper is two fold: to propose a generalized rate model for 
neuron ensembles and to study its property with the use of FPE and the AMM. The paper 
is organized as follows. In Sec. 2, we discuss the generalized rate model for a single cluster 
containing N neurons, investigating its stationary and dynamical properties. In Sec. 3, 
our rate model is extended and applied to an ensemble containing multiple neuron clusters. 
In particular, we study the two-cluster ensemble consisting of excitatory and inhibitory 
clusters. The final Sec. 4 is devoted to discussion and conclusion. 
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2 Single neuron clusters 
2.1 Generalized rate model 

We have adopted a neuronal cluster consisting of N neurons. The dynamics of firing rate 
f i 0) of a neuron i is assumed to be described by the Langcvin model given by 

= F{n)+H{u^)+aG{riMt)+/3^,{t), {i ^ 1 - N) (10) 

with 

= (I) E rkit)+hit), (11) 

where F{x), G{x) and H{x) are arbitrary functions of x: Z {= N ~ 1) denotes the coordi- 
nation number: w is the coupling strength: Ii{t) expresses an input from external sources: 
a and /3 are the strengths of additive and multiplicative noises, respectively, given by rji(t) 
and ^i(t) expressing zero-mean Gaussian white noises with correlations given by 

<T^,{t)rj,{t')> = 6,jS{t-t'), (12) 
<m^At')> = 5,,5{t-t'), (13) 
<V^{t)iAt')> = 0. (14) 

The rate models proposed so far have employed F{x) — —Ax and G{x) — (no multi- 
plicative noises). In this paper, we will adopt several functional forms for F{x) and G{x). As 
for the gain function H{x), two types of expressions have been adopted. In the first category, 
the sigmoid function such as H{x) = tanh(a;), 1/(1 -l-e""^), atan(x), etc. have been adopted. 
In the second category, H{x) is given by the f — I function as H{x) — (x — Xc)Q{x — Xc) 
which expresses the frequency / of autonomous oscillation of a neuron for the applied dc 
current /, Xc denoting the critical value and Q{x) the Heaviside function: 0(a;) = 1 for 
X > and otherwise. It has been theoretically shown in Ref. |65j that when spike inputs 
with the mean ISI (Ti) are applied to an HH neuron, the mean ISI of output signals (To) 
is To — Ti for Ti ^ 15 ms and To ~ 15 for Ti > 15 ms. This is consistent with the recent 
calculation for HH neuron multilayers, which shows a nearly linear relationship between the 
input (ri) and output rates (ro) for n < 60 Hz [36]. It is interesting that the — relation 
is continuous despite the fact that the HH neuron has the discontinuous first-type f — I 
relation. We will adopt, in this paper, a simple expression given by [66j 

Hix) -F^. (15) 

Vx"^ + 1 

although our result is valid for an arbitrary form for H(x). The nonlinear, saturating 
behavior in H{x) arises from the property of the refractory period (r^) where spike outputs 
are prevented for tf<t<tf + Tr after firing at t = tf. 
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2.2 Stationary property 
2.2.1 Distribution of r 

The Fokker-Planck equation for the distribution of p{{ri}, t) is given by [67j 

„2 



drk 
& 



2 V ^"-k 



(16) 



where G'{x) = dG{x)/dx, and 0=1 and in the Stratonovich and Ito representations 
respectively. 

The stationary distribution p{r) for w = and Ii{t) = / is given by 



Inp(r) oc X{r) + Y{r) - ( 1 - | ) In 



with 



X{r) = 2 J dr 
Y{r) 



Fir) 



2 / dr 



a2G(r)2 + /32 
H{I) 



_a^G{r)^+(3^ 

Hereafter we mainly adopt the Stratonovich representation. 
Case I F{x) = —Xx and G{x) = x 

For the linear Langevin model, we get 



ith 



p(r) oc 



Y{r) 



1 



f2H\ 






arctan ^"^^ 


\al3 J 





where H — H{I). In the case oi H = Y{r) = 0, we get the g-Gaussian given by 

p{r) cx [1 - (1 - q)-fr^] ^ , 

with 

7 = 

q = 



2A + a2 
2A + 3a2 



2A + a2 ■ 

We examine the some limiting cases of Eq. (20) as follows. 

(A) For a = and /? 7^ 0, Eq. (20) becomes 

p[r) (X e 1^ . 

(B) For /3 = and a ^ 0, Eq. (20) becomes 

p(r) oc ^-(2A/«^ + l)g-(2H/a^)A^ 



(17) 

(18) 
(19) 



(20) 

(21) 

(22) 

(23) 
(24) 



(25) 



(26) 
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Distributions p{r) calculated with the use of Eqs. (22)-(26) are plotted in Figs. 2(a)- 
2(c). The distribution p{r) for a — 0.0 (without multiplicative noises) in Fig. 2(a) shows 
the Gaussian distribution which is shifted by an applied input / = 0.1. When multiplicative 
noises are added (a 7^ 0), the form oi p{r) is changed to the q-Gaussian given by Eq. (22). 
Figure 2(b) shows that when the magnitude of additive noises (3 is increased, the width 
of p(r) is increased. Figure 2(c) shows that when the magnitude of external input / is 
increased, p{r) is much shifted and widely spread. Note that for a — 0.0 (no multiplicative 
noises), p{r) is simply shifted without a change in its shape when increasing /. 
Case II F{x) = -Ace" and G{x) ^ x'' (a, 6 > 0) 

The special case of a = 1 and = 1 has been discussed in the preceding case I [Eqs. 
(22)-(26)]. For arbitrary a > and 6 > 0, the probabihty distribution p{r) given from Eqs. 
(17)-(19) becomes 

1-1/2 



p{r) (X 



1 



„2b 



exp[X(r)+y(r)], 



with 



X(r) 
Y{r) 



2Ar'^+i 
/32(a + l) 



^ , a + 1 a + 1 
F\ l,-r^,-r^ + l; 



f2Hr\ 



1 1 
2&' 26 



26 
1;^ 



' 26 



(27) 

(28) 
(29) 



where -F(a, 6, c; z) is the hypergeometric function. Some limiting cases of Eqs. (27)- (29) are 
shown in the following. 

(a) The case oi H = Y{r) = was previously studied in Ref. [44] . 

(b) For a = and (3 ^ 0, we get 



p{r) cx exp 



2A 



cx r ^'^/^ exp 



2Hr 



2H\ 



/32 



for a + 1 7^ 
for a -I- 1 = 



(30) 
(31) 



(c) For (i — Q and a 7^ 0, we get 



p{r) (X r ^ exp 



2A 



a2(a-26+l) 



a-26+l 



2H 



,,-2b+l 



^a2(26- 1)^ 
for a - 26 + 1 7^ 0, 26 - 1 7^ (32) 



CX r(2«/" -i/2)exp 



2H 



a2(26- 1) 
2A 



„-2fc+l 



for a - 26 + 1 = (33) 
for 26- 1 = (6 = 1/2) (34) 



cx r 



-[2(\-H)/a^ + l/2\ 



for a - 26 + 1 = 0,26 - 1 = (a = 0, 6 = 1/2) 



(35) 



(d) In the case of a = 1 and 6 = 1/2, we get 



p(r) cx r H 



2w {2\f3^/a*+2H/a^-l/2) 



exp 



(36) 
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which reduces, in the Umit of a = 0, to 



p[r) (X exp 



H 



for a = 



(37) 



Case III F{x) 
We get 



^Alnx and G{x) = x^/^ > 0) 



p{r) oc 



-1/2 



exp 







(^) 









for /? = 



(38) 



Figure 3(a) shows distributions p{r) for various a with fixed values of / = 0.1, b = 1.0, 
a = 1.0 and f3 = 0.0 (multiphcative noise only). With decreasing a, a peak of p{r) at 
r ~ 0.1 becomes sharper. Figure 4(a) shows distributions p{r) for various b with fixed 
values of / = 0.1, a = 1.0, a = 1.0 and /3 = 0.0 (multiplicative only). We note that a change 
in the 6 value yields considerable changes in shapes of p{r). Figures 3(b) and 4(b) will be 
discussed shortly. 

2.2.2 Distribution of T 

When the temporal ISI T is simply defined by T = 1/r, its distribution 7r(T) is given by 



1\ 1 



TTiT)=pl^- j ^. 
For F{x) = -Xx, G{x) = x and /3 = 0, Eq. (26) or (33) yields 



7r(r) cx r(2V"^-i) 



exp 



2H 

„,2 



(39) 



(40) 



which expresses the gamma distribution [see Eq. (7)] [42l|49]. For F{x) — —Xx^, G{x) — 
and /3 = 0, Eq. (32) yields 

7r(r) cx T^^ exp 



2H 



T 



2X\ l_ 
a^) T 



which is similar to the inverse Gaussian distribution [see Eq. (8)] [50]. For F{x) 
G{x) = and /3 = 0, Eq. (38) leads to 



7r(r) cx T-3/2gxp 



(41) 
Inx, 

(42) 



which is similar to the log- normal distribution [see Eq. (9)] |51| . 

Figures 3(b) and 4(b) show distributions of T, tt(T), which are obtained fromp(r) shown 
in Figs. 3(a) and 4(a), respectively, by a change of variable [Eq. (39)]. Figure 3(b) shows 
that with increasing a, the peak of 7r(T) becomes sharper and moves left. We note in Fig. 
4(b) that the form of tt{T) significantly varied by changing b in G{x) — x^ . 

2.2.3 Distribution of R 

When we consider global variables R{t) defined by 



(43) 
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the distribution P{R,t) for R is given by 

j 

Analytic expressions of P{R) are obtainable only for limited cases, 
(a) For f3 ^ and a = 0, P{R) is given by 



P{R) oc exp 



AiV 



where H = H{I). 

(b) For = 0, we get [62] 



with 



1 f 
27r 



dk e'""'' $(fc), 



$(fc) 



where (/>(A:) is the characteristic function for p{r) given by [68] 



m 



-ikx 



p{x)dx, 



21-. (AMfcir I k I), 



with 



A' 



1 

a 



(44) 



(45) 



(46) 
(47) 

(48) 
(49) 

(50) 
(51) 



K^(x) expressing the modified Bessel function. 

Some numerical examples of P{R) are plotted in Figs. 5, 6 and 7 [52). Figures 5(a) and 
5(b) show P{R) for a = 0.0 and a = 0.5, respectively, when N is changed. For a = 0.0, 
P(i?) is the Gaussian whose width is narrowed by a factor of l/\/]V with increasing N. 
In contrast, P{R) for a — 0.5 is the non-Gaussian, whose shape seems to approach the 
Gaussian as increasing TV. These are consistent with the central-limit theorem. 

Effects of an external input / on p{r) and P{R) are examined in Figs. 6(a) and 6(b). 
Figure 6(a) shows that in the case of a = 0.0 (additive noise only), p{r) and P{R) are simply 
shifted by a change in /. This is not the case for a ^ 0.0, for which p{r) and P{R) are 
shifted and widen with increasing /, as shown in Fig. 6(b). 

Figures 7(a) and 7(b) show effects of the coupling w on p{r) and P{R)- For a — 0.0, p{r) 
and P{R) are little changed with increasing w. On the contrary, for a = 0.5, an introduction 
of the coupling significantly modifies p{r) and P{R) as shown in Fig. 7(b). 
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2.3 Dynamical property 
2.3.1 AMM 

Next we will discuss the dynamical property of the rate model by using the AMM. Moments 
of local variables are defined by 



(rf) = 



n,dnpi{n},t)rt (fc-1,2,-) (52) 



Equations of motions of means, variances and covariances of local variables (r^) are given 
by 

'^^'■^^ - (F(rO) + (F(u,)) + ^(G'(r.)G(rO), (53) 



din Tj 



dt N V W, N V 2 

= (r, F{r,)) + (r, Fin)) + {n H{u,)) + (r, H{u,)) 



dt 



2 .nG\r,)G{r,))^'^{r,G'{n)G{u)) 
+ [a^(G(r,)'> (54) 

Equations of motions of the mean, variance and covariance of global variables (i?) are 
obtainable by using Eqs. (43), (53) and (54): 

d{R) _ 1 d{n) 

- iV^^' ^^^^ 

i 

d{B?) 1 ^ ^ d{r, r 



dt ^ ^ dt 

In the AMM [5^, we define fj,, 7 and p given by 



i 

p = {{R-fif), (59) 



^ ^ N 



where fi expresses the mean, 7 the averaged fluctuations in local variables (r^) and p fluc- 
tuations in global variable (R). Expanding in Eqs. (53)- (56) around the average value of 
p, as 

ri = p + 6r^, (60) 

and retaining up to the order of < SriSrj >, we get equations of motions for p, 7 and p 
given by 

^ = /o + /27 + ho + [50^1 + 3(5152 + 5053)7], (61) 

d'y „ „ „, / wN\ ( 7 

it = ^^^^+^'^^(— 

+ (0+l)(.g2 + 25o52)a'7 + a'ffo+/3', (62) 
^ = 2/ip + 2/ii^«p+(</) + l)(.g2 + 23032) «'p+^(a2g2+/32)^ (63) 
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where 

u = wfi + I. (67) 

Original A''-dimensional stochastic equations given by Eqs. (10), (11) and (15) are trans- 
formed to the three-dimensional deterministic equations given by Eqs. (61)-(63). 

Before discussing the dynamical property, we study the stationary property of Eqs. (61)- 
(63). In order to make numerical calculations, we have adopted 

F{x) = -Xx, (68) 
G{x) = X, (69) 

where A stands for the relaxation ratio. Equations (61)-(63) are expressed in the Stratonovich 
representation by 

da , , a^a 

it = -A/^ + /^o + ^, (70) 

^ = _2A, + ^(p-^)+2«^ + «V+/?^ (71) 

I = -2Xp + 2h,wp + 2a'p+'^ + ^^, (72) 

where ho = uj^/v? + 1, hx = l/{u^ + if/^ and h2 = -(3u/2)/(u^ + if'^- The stability of 
the stationary solution given by Eqs. (70)-(72) may examined by calculating eigenvalues of 
their Jacobian matrix, although actual calculations are tedious. 

Figure 8 shows the N dependences of 7 and p in the stationary state for four sets of 
parameters: {a,(i,w) = (0.0,0.1,0.0) (solid curves), (0.5, 0.1, 0.0) (dashed curves), (0.0, 
0.1, 0.5) (chain curves) and (0.5, 0.1, 0.5) (double-chain curves), with /? = 0.1, A = 1.0 and 
/ = 0.1. We note that for all the cases, p is proportional to N~^, which is easily realized in 
Eq. (72). In contrast, 7 shows a weak A'' dependence for a small A'' (< 10). 

2.3.2 Response to pulse inputs 

We have studied the dynamical property of the rate model, by applying a pulse input given 

by 

i{t) = Ae{t-ti)eit2-t) + i^^\ (73) 

with A = 0.5, h = 40, t2 = 50 and = 0.1 expressing the background input, where Q{x) 
denotes the Heaviside function: G(a;) = 1 for a; > and otherwise. 

Figures 9(a), 9(b) and 9(c) show the time dependence of p, 7 and p when the input 
pulse I{t) given by Eq. (73) is applied: solid and dashed curves show the results of AMM 
and DS averaged over 1000 trials, respectively, with a = 0.5, jS = 1.0, A'' = 10 and w = 0.5 
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[69j . Figures 9(b) and 9(c) show that an apphed input pulse induces changes in 7 and p. 
This may be understood from 2a^ terms in Eqs. (71) and (72). The results of AMM shown 
by solid curves in Figs. 9(a)-(c) are in good agreement with DS results shown by dashed 
curves. Figure 9(d) will be discussed in the foUowings. 

It is possible to discuss the synchrony in a neuronal cluster with the use of 7 and p 
defined by Eqs. (58) and (59) 56J. In order to quantitatively discuss the synchronization, 
we first consider the quantity given by 

= E < [^^(*) - ^^■(^)]' 2[7(0 - p{t)]. (74) 



When all neurons are in the completely synchronous state, we get ri{t) = R{t) for all i, and 
then P{t) = in Eq. (74). On the contrary, we get P{t) ^ 2(1 - 1/A^)7 = Po(^) in the 
asynchronous state where p ~ 7/iV [56j . We may define the synchronization ratio given by 

which is and 1 for completely asynchronous (P = Pq) and synchronous states (P — 0), 
respectively. Figure 9(d) shows the synchronization ratio S{t) for '-f{t) and p{t) plotted 
in Figs. 9(b) and 9(c), respectively, with a ~ 0.5, (3 — 1.0, = 10 and w = 0.5. The 
synchronization at t < 40 and t > 60 is 0.15, but it is decreased to 0.03 at 40 < t < 50 
by an applied pulse. This is because 7 is more increased than p by an applied pulse. The 
synchronization ratio is vanishes for w = 0, and it is increased with increasing the coupling 
strength 56J. 

Next we show some results when indices a and b in F{x) = —Ax" and G{x) — are 
changed. Figure 10(a) shows the time dependence of p for (a, 6) = (1, 1) (solid curve) and 
{a,h) = (2,1) (dashed curve) with a = 0.0, /3 = 0.1, iV = 10 and w = 0.0. The saturated 
magnitude of p for a — 0.5 is larger than that for a = 0.0. Solid and dashed curves in 
Fig. 10(b) show p for (a, 6) = (1,1) and (1,0.5), respectively, with a = 0.5, /3 = 0.001, 

= 10 and w = 0.0. Both results show similar responses to an applied pulse although 
p for a background input of 7^''' = 0.1 for (a, 5) — (1,0.5) is a little larger than that for 
(a,fe) = (l,l). 

2.3.3 Response to sinusoidal inputs 

We have applied also a sinusoidal input given by 



m - A 



/27rA 



(76) 



with A = 0.5, = 0.1, and Tp = 10 and 20. Time dependences of p for Tp = 20 and 
Tp = 10 are plotted in Figs. 11(a) and 11(b), respectively, with a — 0.5, [3 — 1.0, w ~ 0.0 
and A = 10, solid and dashed curves denoting p and /, respectively. The delay time of p 
against an input I(t) is about tj^ ^ 1.0 independent of Tp. The magnitude of p for Tp — 10 
is smaller than that for Tp = 20. 
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3 Multiple neuron clusters 
3.1 AMM 

We have assumed a neuronal ensemble consisting of M clusters, whose mth cluster includes 
A^„j neurons. The dynamics of firing rate (> 0) of a neuron i in the cluster m is assumed 
to be described by the Langevin model given by 

^^^^ = F{rmi) + H{Umi) + amG{rrni)rirni{t) + Pmivai{t), (77) 

(m = 1 - A/, i = 1 - Nra) 

with 

+ Im{t), (78) 

where F(x) and G{x) are arbitrary functions of x: H{x) is given by Eq. (15): Z.„-i (= A'^,,,, — 1) 
denotes the coordination number: Im{t) expresses an external input to the cluster m: am 
and (3m are the strengths of additive and multiplicative noises, respectively, in the cluster 
m given by 7]mi (t) and ^mi (t) expressing zero- mean Gaussian white noises with correlations 
given by 

<Vm{t)l]mj{t') > SnmS^JS{t~t'), (79) 

< Cm{t) Uj{t') > = 5nm5rAi-t')^ (80) 

< Vm{t) > = 0. (81) 

In the AMM [56j . we define means, variances and covariances, /x™, 7^ and Pmm given 

by 



I 
z 

Pmn = ((^m - Mm)(-Rn - Mn)), (84) 

where the global variable Rm is defined by 

I 

Details of deriving equations of motions for /i„i, 7^ and are given in the Appendix 
[Eqs. (A10)-(A12)]. 

3.2 Stationary property 

Now we consider an E-I ensemble (M = 2) consisting of excitatory (E) and inhibitory (I) 
neuron clusters, for which we get equations of motions for jm and pmn from Eqs. 
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(A10)-(A16): 



Ie.O + fE,2lE + h 



E,0 



[gE,o9E,i + 3{gE,igE.,2 + gE,ogE.,3)iE], 



dm 

dt 



—TT = //,0 + fiXII + ^Ifi 



[gi.ogiA + 3(57437,2 + gi,ogi,3hi] 



(86) 



(87) 



dt 



2/7;, 17b + 2/17;,! 



Y weeNe\ f IE 
i.4> + l)(5l.i + 2gEfigE.2)aElE + alsl.o + 



weiPei 



^ - 2/7,177 + 2/17,1 



-wiiNi \ 
Zi J 



WiePEI 



dp 



EE 



dt 



dpii 

dt 



(0 + l)igl^ + 2g7,o57,2)a?7/ + a?5?,o + Pi 

'2fE,lPEE + 2hE,liwEEPEE ~ WeiPei) 

("Isio + Pi) 



(-^ + l)(57J,i + '2gE,ogE,2) aE Pee + 



^fi,iPii + 2hi^i{-~wiipii + wiePei) 
(</* + l)(57,i + 257,o5/,2) aj pu + 



(ajglo 



PI 



dp 



EI 



dt 



{fE,l + fl,l)PEI 
hE,l{wEEPEI - WElPll) 



Ni 



hi,i{-wiipEi + wiePee) 



(89) 



(90) 



(91) 



1) 



[(.9l,i + 2.g7;,o5B,2) a% + (fff.i + 257,057,2) a]]pEi- 



(92) 



Here we set — W77 < and —wei < after convention. Equations (86)-(92) correspond to a 
generalized Wilson-Cowan (WC) model, because they reduce to WC model HH] if we adopt 
F{x) = —A and G{x) = x (see below), neglecting all fluctuations of 7^ and p^^' (77, if^E, I). 

It is difficult to obtain the stability condition for the stationary solution of Eqs. (86)- 
(92) because they are seven-dimensional nonlinear equations. We find that equations of 
motions oi he and p,i are decoupled from the rest of variables in the cases of G{x) — x and 
G{x) = x^l"^ with F(x) = -A x, for which /^,2 = and (5r,,i5r/,2 + gr,,og7j,3) = iv = 
in Eqs. (86) and (87). We will the discuss the stationary solutions for these two cases in 
the followings. 

(A) F{x) = -Aa; and G{x) = x 

Equations of motions for pE and pi become 

dpE 



dt 

dpi 
dt 



= - Xe- 



Xr 



PE + H{weEPE - WEI pi + Ie): 



Pi + H{wiEPE - ^77^/ + Ii)- 



(93) 
(94) 
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The stationary slution is given by 

= fl{llE, ^.l) = - [\e fJ-E + H{wEEtJ'E ~ WEIfJ.1 + Ie), (95) 

= f2{fJ.E, f^i) = - - ^J'I + H{wiEf^E - wiim + Ii). (96) 
The stabiHty condition for stationary solutions is given by 

2 2 

T = -Ab - A/ + ^ + ^ - WEEhEA - wiihi^i < 0, (97) 

D = (^Xe ^-^^ WEEhE.ij ^A/ - - wiihis 

+ WEiwiEhE,ihiA > 0, (98) 

where T and D denote the trace and determinant, respectively, of Jacobian matrix of Eqs. 
(93) and (94). 

By solving Eqs. (93) and (94), we get stationary solutions of ^e and /i/. Figures 12(a) 
and 12(b) show ^e said fij, respectively, as a function of wee for various a {= oe — cti) 
with Xe = Xj = 1.0, Wei = wje — wu = 1, = // = and Ne ^ Nj = 10. In the case of 
a = 0, fiE and /i/ are zero for wee < Wc but become finite for wee > Wc where Wc (=1.5) 
denotes the critical couplings for the ordered state. With increasing a, the critical value of 
Wc is decreased and magnitudes of /i^ and /i/ in the ordered state are increased. This shows 
that multiplicative noise works to create the ordered state |39j . 

(B) F{x) ^ -Xx and G{x) = x^/^ 

Equations of motions for fiE and /i/ become 



= -XE^iE + ^ + H{WEE^^E - WEIfJ^I + Ie), (99) 

at 4 

+ -f + H{wiE^lE - Wiifii + Ii). (100) 
at 4 

The stationary solution is given by 

Q,2 

= fliHE,tl'l) = -XeUE + ^ + H{wEEtJ-E - WEI^iI + Ie), (101) 
= f2{t^E,fJ'l) ^ -Xifij + ^ + H{wiEt^E - WiJlJ.1 + Ij). (102) 

The stability condition for stationary solution is given by 

T = -Xe - A/ - WEEhE.i - wiih^i < 0, (103) 
D = [Xe - WEEhE,i){Xi - wiihj^i) + WEiwiEhE,ihiA > 0, (104) 

where T and D express the trace and determinant, repectively, of Jacobian matrix of Eqs. 
(99) and (100). 

Figure 13(a) and 13(b) show the wee dependences of fiE and fij, respectively, for various 
a {— aE ~ ai) with A^; = A/ = 1.0, wei = wje = wu ~ 1.0 and Ne = Nj = 10. Equations 
(99) and (100) show that multiplicative noise play a role of inputs, yielding finite ^e and 
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Hi even for no external inputs {Ie — Ii = 0). With increasing a, magnitudes of /x^ and /ii 
are increased and the critical value of Wc for the ordered state is decreased. It is interesting 
to note that the behavior of fiE and fij at wee ^ Wc is gradually changed with increasing 
a. 

3.3 Dynamical property 

We have studied the dynamical property of the rate model, by applying pulse inputs given 
by 

Ir,it) = Ar,e{t-h)eit2-t)+lj^''\ iV = E,I) (105) 

with Ae = 0.5, Aj = 0.3, 4''^ =0.1, /f = 0.05, ti = 40 and = 50. The time courses 
of iirj, "fri and prirj' {VtV' — E , I) arc plotted in Figs. 14(a)-14(c), respectively, where solid, 
chain and dashed curves show the results of AMM and dotted curves those of DS averaged 
over 1000 trials for ue = cti = 0.5, fiE = = 0.1, and wee = wei = wie = wu = 1.0. 
The results of AMM are in good agreement with DS results. 

Responses of he and /x/ for various sets of couplings are plotted in Figs. 15(a)-15(f): 
twlOOl in Fig. 15(b), for example, means that {wee,wei,wie,wii) =(1,0,0,1). Figure 
15(a) shows the result of no couphngs {wee = wei = wje = wn = 0). When intracluster 
couplings of Wee = 1-0 and wu = 1.0 are introduced, /i/ in the inhibitory cluster is much 
suppressed, while the excitatory cluster is in the ordered state with he = 0.73 for no input 
pulses, as shown in Fig. 15(b). Figure 15(c) shows that when only the interclustcr coupling 
of Wei is introduced, the magnitude of he is decreased compared to that in Fig. 15(a). In 
contrast. Fig. 15(d) shows that an addition of only wie enhances magnitude of Hi compared 
to that in Fig. 15(a). We note in Fig. 15(e) that when both intercluster couplings of wei 
and Wie are included, magnitude of /x^ is considerably reduced while that of hi is slightly 
increased. When all couplings are added, magnitudes of both he and hi are increased 
compared to those for no couplings shown in Fig. 15(a). Figure 15(a)-15(f) clearly shows 
that responses of he and hi to inputs significantly depend on the couplings. 

Figure 16(a) and 16(b) show the synchronization ratios of Se and Si, respectively, 
defined by [see Eq. (75)] 

S,it)=[ ''^^^f^2f-' ), irj = E,I) (106) 

for various sets of couplings when inputs given by Eq. (105) are applied (Ae ~ 0.5, Ai = 0.3, 
= 0.1, = 0.05, aE = ai = 0.5, /Se = Pi = 0.1, and Ne = Ni = 10). First we discuss 
the synchronization ratio at the period of t < 40 and t > 60 where the pulse input is not 
relevant. With no intra- and inter-cluster couplings {wee = wei = wie = wn = 0), we 
get Se = Si = 0. When only the intra-cluster couplings of wee = 1 and wn = 1 are 
introduced, we get Se = 0.15 and Si = —0.67, as shown by dashed curve in Fig. 16(b). 
When inter-cluster coupling of Wei = 1 is included, the synchronization in the excitatory 
cluster is decreased to Se = 0.08 (dotted curves in Fig. 16(a). In contrast, when inter-cluster 
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of wiE = 1 is introduced, the synchrony in the inhibitory cluster is increased to Sj = 0.06, 
as shown by dotted curve in Fig. 16(b). When inter-cluster couplings of Wei = Wie = 1 are 
included, the synchronization ratios almost vanish (chain curves). Solid curves show that 
when both intra- and inter-cluster couplings are included, we get Se = 0.24 and Si = 0.04. 
It is noted that the responses of the synchronizations to a pulse applied at 40 < f < 50 are 
rather complicated. When an input pulse is applied at t = 40, the synchronization ratios 
are generally decreased while Se with wOllO increased: 5'^; with wOlOO is once decreased 
and then increased. When an applied pulse disappears at t = 50, the synchronization are 
increased in the refractory period though the synchronization ratios for wOlOO and wOllO 
are decreased. 

Figures 17(a)-17(e) show responses of an E-I ensemble with Ne = Nj = 10 when the 
input pulse shown in 17(a) is applied only to the excitatory cluster. Responses of the local 
rate of of single neurons in the excitatory [r] = E) and inhibitory clusters [r] = I) are 
shown by solid and dashed curves, respectively, in Fig. 17(b). Local rates in Fig. 17(b) 
which are obtained by DS with a single trial, have much irregularity induced by additive 
and multiplicative noises with cue = aj = 0.1 and /3e = /?/ = 0.1. Figure 17(c) shows 
the population-averaged rates of Rri{t) obtained by DS with a single trial. The irregularity 
shown in Fig. 17(c) is reduced compared to that of r^(t) in Fig. 17(b), which demonstrates 
the advantage of the population code. When DS is repeated and the global variable is 
averaged over 100 trials, we get the result of Rr,{t) whose irregularity is much reduced by 
the average over trials, as shown in Fig. 17(d). The AMM results of /ir,(f) shown in Fig. 
17(e) are in good agreement with those shown in Fig. 17(d). 

4 Discussion and conclusion 

We may calculate the stationary distributions in the E-I ensemble. Figures 18(a)-18(c) show 
global distributions of Pe{R) and Pi{R) which are averaged within the excitatory and in- 
hibitory clusters, respectively. We have studied how distributions are varied when couplings 
are changed: «;0110 in Fig. 18(b), for example, means {wee,Wei,wie,wii) = (0,1,1,0). 
Figure 18(a) shows the results without couplings, for which distributions of Pe{R) and 
Pi{R) have peaks at i? ~ 0.1 and 0.05, respectively, because of applied background inputs 
(J^j''-' = 0.1 and Ij''^ = 0.05). When intcrcluster couplings of wei = 1-0 and wje — 1-0 are 
introduced, the peak of Pe{R) locates at i? ~ 0.02, while that of Pi{R) is at i? ~ 0.08: 
their relative positions are interchanged compared to the case of no couplings shown in Fig. 
18(a). When all couplings with wee =wei = wje = wu = 1.0 are included, we get the dis- 
tributions shown in Fig. 18(c), where both Pe{R) and Pi{R) have wider distributions than 
those with no couplings shown in Fig. 18(a). Our calculations show that the distributions 
are much influenced by the magnitudes of couplings. 

We have proposed the rate model given by Eqs. (10) and (11), in which the relaxation 
process is given by a single F{x). Instead, when the relaxation process consists of two terms: 

F{x) ^ ciFi{x) + C2F2{x), (107) 
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with ci + C2 — 1, the distribution becomes 



p{r)^[p,{r)r[p,{r)p, (108) 

where Pk{r) {k — 1,2) denotes the distribution only with F{x) = Fi{x) or F{x) = Fi{x). 
In contrast, when multiplicative noises arise from two independent origins: 

axri{t) — > ciaixr]i{t) + 0202x772 (i), (109) 

the distribution for /3 = H = becomes 

p(0«r~''^^^''"'+''"^^+''- (110) 
Similarly, when additive noises arise from two independent origins: 

m)^Cif3iUt)+C2p2Ut), (111) 

the distribution for a = H = becomes 

p(r) cxe-^/("i"i+"^"^). (112) 

Equations (108), (110) and (112) are quite different from the form given by 

p{r) = cipi{r) + C2P2{r), (113) 

which has been conventionally adopted for a fitting of theoretical distributions to that 
obtained by experiments. 

It is an interesting subject to decode the stimulus from the observed spiking rate of 
neurons. In BMI, stimulus y{t) is decoded from observed rate signals r{t — u) with the use 
of Eq. (4). In the AMM, the relation between the input I{t) and the averaged rate of R{t), 
fi{t), given by Eq. (70) yields that I{t) is expressed in terms of fi{t) and d^{t)/dt as 

m = [^^M+(A^£ad^-MO, (114) 

^l-[d^i/dt+{\-a-^/2)ii{t)f 
~ ^ + ( A- ^ - w ) ^(i). forsman/x(0 (115) 



dt \ 2 

The dfj,/dt term plays an important role in decoding dynamics of iJ,(t). In an approach using 
the Bayesian statistics |70) . the conditional probability of an input / for a given rate i?, 
P{I I i?), is expressed by 

P{I \ R) (X P{R\ I) P{I), (116) 

where P{R \ I) is the conditional probability of R for a given / and P{I) the likelihood of 
/. An estimation of the value of / which yields the maximum of P{I \ R), is known as the 
maximum a posteriori (MAP) estimate. Since P{R \ I) is obtainable with the use of our 
rate model, as shown in Figs. 6(a) and 6(b), we may estimate P{I \ R) by Eq. (116) if 
P{I) is provided. We note that I{t) given by Eq. (114) corresponds to the center of gravity 
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of P{I I R), which is expected to be nearly the same as the maximum value obtained 
by MAP estimate. More sophisticate Bayesian approach using the recursive method has 
been proposed [7D], although it is not unclear whether such Bayesian networks may be 
implemented by real neurons. 

The structures of neuronal networks have been discussed based on the theory on complex 
networks [71] [72] . The neural network of nematode worm C. elegans is reported to be small- 
world network. It has been recently observed by the functional magnetic resonance imaging 
(fMRI) that the functional connectivity in human brain has the behavior of the scale-free 
[73] or small- world network [74] [75]. Most of theoretical studies have assumed the local 
or all-to-all coupling in neuron networks, as we have made in our study. In real neural 
networks, however, the couplings are neither local nor global. A new approach extending 
the AMM has been proposed to take into account the couplings from local to global and/or 
from regular to random couplings [60] . It has been shown that the synchronization in the 
small-world networks is worse than that in the regular networks due to the randomness 
introduced in the small-world networks. 

In recent years, it becomes increasingly popular to study the distributed information 
processing by using cultured neuronal networks which are cultivated in an artificial way 
[76j . It is possible that every cell in the cultured network may be observed, monitored, 
stimulated and manipulated with high temporal and spatial resolutions. The observed ISI 
distributions of the cultured networks with 50-10^ neurons are reported to obey the scale-free 
distribution 77J. Although our AMM study in Sec. 3 has been made for an E-I (M = 2) 
cluster, it would be interesting to investigate the dynamics of larger ensembles modeling 
complex networks and cultured networks, which is left for our future study. 

To summarize, we have discussed the stationary and dynamical properties of the general- 
ized rate model by using the FPE and AMM. The proposed rate model is a phenomenological 
one and has no biological basis. Nevertheless, the generalized rate model is expected to be 
useful in discussing various properties of neuronal ensembles. Indeed, the proposed rate 
model has an interesting property, yielding various types of stationary non-Gaussian dis- 
tributions such as gamma, inverse-Gaussian and log-normal distributions, which have been 
experimentally observed [47j-|51j. The stationary distribution and dynamical responses of 
neuronal clusters have been shown to considerably depend on the model parameters such 
as strengths of noises and couplings. A disadvantage of our AMM is that its applicability 
is limited to weak-noise cases. On the contrary, an advantage of the AMM is that we can 
easily discuss dynamical property of an iV-unit neuronal cluster. In DS and FPE, we have 
to solve the A^-dimensional stochastic Langevin equations and the {N + l)-dimensional par- 
tial differential equations, respectively, which are more laborious than the three-dimensional 
ordinary differential equations in the AMM. We hope that the proposed rate model in the 
AMM is adopted for a wide class of study on neuronal ensembles. 
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APPENDIX 

A AMM for multiple clusters 

We will present a detail of an application of the AMM to the rate model describing multiple 
clusters given by Eqs. (77) and (78). The Fokker-Planck equation for the distribution of 
Pilrmi}, t) is given by [67] 
i9 , , 

= -Yl 7r^{[^(^™'=) + ^G'(r„fe)G(r™,) + Hiu^k)] p({r™}, t)} 

mk 



2'^ d' 

■mk 



ik 



where G'{x) = dG{x)/dx, and 0=1 and in the Stratonovich and Ito representations, 
respectively. 

When we consider global variables of the cluster m given by 

the distribution P(Rm,t) for i?„j is given by 

P{Rra,t) = /"•••/"n»dr™p({r™},t)5(i?„- ^^r„j). (A3) 

Variances and covariances of local variables are defined by 

{'rt^'rn,) = / n„u dr^, p({r„„}, t) rf„, r^', . (fc, k' = 1, 2, ••) (A4) 



Equations of motions of means, variances and covariances of local variables (r^i) are given 

by 



dt 



= (i^(r™)) + (^("m.)> + V^(G'(r„„)G(r„,)), (A5) 



+ — ^^{'''miG' {rnj)G{rnj)) H ^^{fnjG' (rmi)G{rmi)) 

(A6) 

Equations of motions of the mean, variance and covariancc of global variables (i?m) are 
obtainable by using Eqs. (A3), (A5) and (A6): 

d{Rm) _ 1 d{r,ni) ... 

dt Nrn ^ dt ' ^ ' 
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d{R„iRn 
dt 



N N 



» 3 



d(j'mi ^nj) 

di ■ 



(A8) 



Variances and covariances, fimi Im and pmn are given by Eqs. (82)- (84). Expanding rj, 
in Eqs. (A5)-(A8) around the average value of jirn ^-s 



(A9) 



and retaining up to the order of < SrmiSrmj >, we get equations of motions for fx^, 7to and 
Pmn given by 

dlJ.m 



dt 



dOm 

dt 



dt 



+ 



fm.O fm., 

'<j)al ■ 



[gmftgm,! + 3(5m,lfl'm,2 + 5m,05m,3)7n 



+ 



M 



(AlO) 



m, 

m.O 



^mnPmn 

(All) 



(/m,l + fn,l)pmn + ^m,l 



M- 1 / ^ 



'^mn' Pnn' 



M-l/ ^ 



'^nn' Pmn* 



n'(^n) 



+ ^^^-^[(fi'm,! + 25m,0fi'm,2) Q;^ + (flll + "^gnfignfi) al]p„ 



nZ. 



(A12) 



where 



fm,e 
gm,e 



1 d^F{nm) 
1 d^H{u^) 



M 



For a two-cluster (M = 2) ensemble consisting of excitatory and inhibitory clusters, equa- 
tions of motions given by Eqs. (A10)-(A12) reduce to Eqs. (86)-(92). 
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Figure 1: Response of a 10-unit FN neuron cluster subjected to additive and multiplicative 
noises; (a) input signal I{t), (b) a local membrane potential [v{t)] and (c) a global membrane 
potential [V{t) = (l/N) Vi{t)] obtained by direct simulation (DS) with a single trial: (d) 
a global membrane potential [V{t)] calculated by DS with 100 trials: (e) the result of the 
AMM \ji{t)]. 

Figure 2: (a) Distributions p{r) of the (local) firing rate r for various a with A = 1.0, 
l3 = 0.1, / = 0.1 and w = 0.0, (b) p{r) for various /3 with A = 1.0, a = 1.0, I = 0.1 and 
w = 0.0, and (c) p{r) for various / with A = 1.0, a = 0.5, /3 = 0.1 and w = 0.0. 

Figure 3: (a) Distributions p{r) of the (local) firing rate r and (b) 7r(T) of the ISI T for 
a = 0.8 (chain curves), a = 1.0 (solid curves), a = 1.5 (dotted curves) and a = 2.0 (dashed 
curves) with b = 1.0, A = 1.0, a = 1.0, /3 = 0.0 and / = 0.1. 

Figure 4: (a) Distributions p{r) of the (local) firing rate r and (b) 7r{T) of the ISI T for 
b = 0.5 (dashed curves), b = 1.0 (solid curves), b = 1.5 (dotted curves) and b = 2.0 (chain 
curves) with a = 1.0, A = 1.0, a = 1.0, (3 = 0.0 and / — 0.1: results for b — 1.5 and 6 = 2 
should be multiplied by factors of 2 and 5, respectively. 

Figure 5: Distributions P{K) of the (global) firing rate R for (a) a = 0.0 and (b) a = 0.5, 
with iV = 1, 10 and 100: A = 1.0, /3 = 0.1, w = 0.0 and I = 0.1. 

Figure 6: Distributions p{r) (dashed curves) and P{R) (solid curves) for (a) a = 0.0 and 
(b) a = 0.5 with I = 0.1 and / = 0.2: iV = 10, A = 1.0, /3 = 0.1 and w = 0.0. 

Figure 7: Distributions p{r) (dashed curves) and P(i?) (solid curves) for (a) a = 0.0 and 
(b) a = 0.5 with w = 0.0 and w = 0.5: A/' = 10, A = 1.0, ^ = 0.1 and / = 0.1. 

Figure 8: The N dependence of 7 and p in the stationary states for four sets of parameters: 
(a, /?,«;) = (0.0,0.1,0.0) (solid curves), (0.5,0.1,0.0) (dashed curves), (0.0,0.1,0.5) (chain 
curves) and (0.5,0.1,0.5) (double-chain curves): A = 1.0, A' = 10 and I = 0.1. 

Figure 9: Time courses of (a) /x(t), (b) 7(t), (c) p(t) and (d) 5(i) for a pulse input /(i) 
given by Eq. (73) with A = 1.0, a = 0.5, (3 = 0.1, A' = 10 and w = 0.5, solid and chain 
curves denoting results of AMM and dashed curves expressing those of DS result with 1000 
trials. 

Figure 10: (a) Response of /i(t) to input pidse I{t) given by Eq. (73) for (a, 6) = (1,1) 
(solid curve) and (a, 6) = (2, 1) (dashed curve) with a = 0.0, /3 = 0.1, A^ = 10 and A = 1.0. 
(b) Response of fi{t) to input pulse I{t) for (a, 6) = (1, 1) (solid curve) and (a, 6) — (1, 0.5) 
(dashed curve) with a = 0.5, = 0.001, AT = 10, A = 1.0 and w = 0.0. 

Figure 11: Response of ^{t) (solid curves) to sinusoidal input I{t) (dashed curves) given 
by Eq. (73) for (a) Tp = 20 and (b) Tp = 10 with A = 0.5, A = 1.0, a = 0.5, /? = 0.1, u; = 
and AT = 10 (a = 1 and 6 = 1). 

Figure 12: Stationary values of (a) he and (b) /xj in an E-I ensemble as a function of 
Wee for various values of a {= aE = cxi) for the case of G{x) = x with A^ = A/ = 1.0, 
Wei = wiE = wii = 1.0 and Ie = Ii = and A'^; = Nj = 10. 
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Figure 13: Stationary values of (a) /ie and (b) /U/ in an E-I ensemble as a fiuietion of 
Wee for various values of a {= ue = cti) for the case of G{x) = x^l"^ with = A/ = 1.0, 
Wei = wiE = wii = 1.0 and Ie = Ii = and Ne = TV/ = 10. 



Figure 14: Time courses of (a) ije and /ij, (b) 7_e and 7/, and (c) pee, Pii and pEi, 
for pulse inputs given by Eq. (105) with cue = ai {= a) = 0.5: (3e = Pi {= P) = 0.1, 
Ne = iV/ = 10 and wee = wei = wie = wn = 1.0: solid, dashed and chain curves denote 
the results of AMM and dotted curves express those of DS with 1000 trials. 



Figure 15: Responses oi he (solid curves) and jii (dashed curves) of an E-I ensemble 

to pulse inputs given by Eq. (105); (a) (t<;_E_E, w^/, w/b, w//) = (0, 0,0,0), (b) (1,0,0,1), (c) 
(0,1,0,0), (d) (0,0,1,0), (e) (0,1,1,0) and (f) (1,1,1,1): wlOOl, for example, expresses the case 

of (b): as = ai = 0.5, (3e = Pi = 0.1, Ae = 0.5, Ai = 0.3, jf^ = 0.1, if^ = 0.05 and 
Ne = Ni = 10. 



Figure 16: Synchronization ratios of (a) Se and (b) Si of an E-I ensemble to pulse inputs 
given by Eq. (105) with Ae = 0.5, Ai = 0.3, ' = 0.1, /f = 0.05, = a/ = 0.5, Pe = 
Pi = 0.1, and Ne = Ni = 10: twlOOl, for example, denotes {wee,wei,wie,wii)={1,0,0,1). 



Figure 17: Responses of an E-I ensemble; (a) input signal lE{t), (b) local rates rrj{t) and 

(c) global rate i?^(t) (rj = E and /) obtained by direct simulation (DS) with a single trial: 

(d) global rates Rr,{t) calculated by DS with 100 trials: (e) the result of the AMM jinit)'- 
ue = ai = 0.5, Pe = Pi = 0.1, wee = wei = wie = wii = 1-0, Ae = 0.5, Ai = 0.0, 

= 0.1, /j^^ = 0.05 and Ne = Nj = 10. Sohd and dashed curves in (b)-(e) denote the 
results for excitatory (E) and inhibitory (I) clusters, respectively. 



Figure 18: Stationary global distributions in E [Pe(B)] and I clusters [Pj(B)] for (a) 
{wEE,WEi,wiE,wii)^{0,Ofl,0), (b) (0,1,1,0) and (c) (1,1,1,1) with I^j^^ = 0.1, /f' = 0.05, 
aE = Oil (= a) = 0.5, Pe = Pi (= P) = 0.1 and Ne = Nj = 10, solid and dashed curves 
denoting Pe{R) and Pi{R), respectively. 
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